The Glass Transition Temperature of Water: A Simulation Study 
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We report a computer simulation study of the glass transition for water. To mimic the difference 
between standard and hyperquenched glass, we generate glassy configurations with different cooling 
rates and calculate the T dependence of the specific heat on heating. The absence of crystallization 
phenomena allows us, for properly annealed samples, to detect in the specific heat the simultaneous 
presence of a weak pre-peak ("shadow transition"), and an intense glass transition peak at higher 
temperature. We discuss the implications for the currently debated value of the glass transition 
temperature of water. We also compare our simulation results with the Tool-Narayanaswamy- 
Moynihan phenomenological model. 

PACS numbers: 



Much recent research has focused on the properties 
of glassy water, the most common form of water in 
the universe, which can exist in more than one distinct 
amorphous formQ],Q,[^. The conversion between differ- 
ent glass structures, the different routes producing glass 
structures, and the relation between the liquid and the 
glass phases are under active debate. 

A particularly relevant aspect of this debate concerns 
the identification of the glass transition temperature T g 
at ambient pressure and the magnitude of the associated 
jump of the specific heat, an issue which has relevance 
also for determining the fragility of water. Extrapolation 
of T g in binary aqueous solutions, in the limit of vanishing 
solute concentration, provides the estimate T g ks 136 K 
0. Early differential scanning calorimetry (DSC) stud- 
ies report conflicting results. Some experiments detect 
the glass transition [j| but others do not . An exother- 
mal peak in the specific heat of properly-annealed h 
perquenched water supports the estimate T g w 136 K 
with a specific heat jump of 1.6 — 1.9 J/mol/K. This 
T g valued, H has been recently debated [nj El 112 • It 
has been suggested [l^ that the small peak measured in 
Ref. is a pre-peak typical of annealed hyperquenched 
samples preceding the true glass transition located at 
Tg w 165 K. Assigning T g w 165 K would explain many 
of the puzzles related to the glass transition in water 
|il Ell [13 Unfortunately, the T g w 165 K proposal can 
not be experimentally tested due to the homogeneous 
nucleation of the crystal phase at T x ks 150 K. 

Here we report a numerical study of the temperature 
dependence of the specific heat across the glass-to-liquid 
transition for the extended simple point charge (SPC/E) 
model for water. We analyze the effects both of the cool- 
ing rate and of annealing ("aging") before heating the 
glas s, since both effects are important for determining T g 
[Tj, [l4( , and both effects have been studied extensively 
in many materials [l5l Il6| . Numerical studies are partic- 



ularly suited since crystallization does not take place on 
the time scale probed in simulations. With an appropri- 
ate choice of the heating and cooling rates to mimic the 
experimental conditions, we show that both the shadow 
and the glass transition peaks can be resolved in the 
same heating scan. Finally, we compare the simulation 
results with the Tool-Narayanaswamy-Moynihan (TNM) 
phenomenological model [l5| uji Eg] ■ 

We perform NVT molecular dynamics (MD) simula- 
tions for a system of N = 216 molecules, with periodic 
boundary conditions. Interactions are cut off at a dis- 
tance of r = 2.5cr, where a is the length parameter 
defined in the SPC/E potential, and the reaction field 
method is implemented to account for the long range in- 
teractions. We average quantities over 32 independent 
trajectories at fixed density p = 1 g/cm 3 . During cool- 
ing or heating, T is continuously changed by ST = qSt, 
where q is the cooling/heating rate, and St — 1 is is 
the elementary time step. We perform: (i) cooling scans 
at constant cooling rate down to T = K, starting 
from equilibrium liquid configurations at T — 300 K, 
(ii) heating scans at constant heating rate, starting from 
T = K glass configurations, (iii) aging at constant 
Tage = 100 K, where significant aging effects are ob- 
served. We study two cooling rates q c = — 3 x 10 10 K/s 
and q c = — 10 13 K/s, to mimic respectively the standard 
and hyperquenched cooling rates, and one heating rate 
qh = +3 x 10 10 K/s. Slow experimental scan rates are 
typically w 0.3 K/s, while the slowest simulation scan 
rate compatible with present computational facilities is 
10 7 times faster («3x 10 10 K/s). Hence, the tempera- 
ture at which the system will lose equilibrium on cooling 
will be significantly higher in simulations than in exper- 
iment. Still, the key fact that the structural relaxation 
time becomes longer than the experimental (or simula- 
tion) time is the same for experiments and simulations. 
Therefore, as we will show below, while the T g estimates 
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differ, the T-dependence and the phenomenology do not 
depend significantly on the scan rate. In hyperquench ex- 
periments, a cooling rate 10 4 times faster than the slow 
or 'standard' rate is usually achieved, while in the present 
simulations the faster quench rate is approximately 300 
times faster than the slower quench rate. 

Figure H shows the specific heat Cy(T) calculated by 
differentiating the temperature dependence of the to- 
tal energy of the system on heating at the rate = 
+3 x 10 10 K/s. The glass configurations are obtained by 
cooling equilibrium T = 300 K liquid configurations at 
the 'standard' cooling rate q c — — 3 x 10 10 K/s. Follow- 
ing the usual experimental protocol, we estimate T g from 
the intersection of the two dashed lines in FigQ The re- 
sulting value, T g = 188K, is slightly below the lowest T 
at which equilibrium simulations can be performed for 
SPC/E[H Hi • The CV(T) rise of « 55 J/mol/K, is 
more than an order of magnitude larger than the exper- 
imentally measured rise of « 1.6 — 1.9 J/mol/K[l3j. For 
T > 240 K, CV(T) coincides with equilibrium data for 
the SPC/E potential [2ll| . Indeed, the equilibrium relax- 
ation time of the system for T > 240 K is < 20 ps, smaller 
than the characteristic scan time 1 K/qh ~ 30 ps. 

We next compare in Figs. Ufa) and|2Jb) the behavior 
of Cy(T) on heating two different glasses, the 'standard 
glass' obtained with the cooling rate q c = — 3 x 10 10 K/s, 
and the 'hyperquenched glass' obtained with the faster 
rate q c = — 10 13 K/s. For the hyperquenched glass, 
Cy(T) develops a valley for T < T„, in agreement with 
DSC heating scan experiments [12. 1221 l23l |24[ Ci ndeed. 
FigEI a ) is remarkably similar to Fig. 1 of Ref . [23] ) . The 
presence of a valley can be related to the descent of the 
system on the potential energy landscape upon heating 
with a rate slower than the cooling rate[2^. Figs. Ufa) 
and 01)) show CV(T) for the heating scan of the hyper- 
quenched glass which has been annealed at T agc = 100K, 
for different aging times < t age < 20 ns. This annealing 
procedure is intended to mimic the experimental anneal- 
ing procedure |25j. Note that aging reduces the valley in 
CV(T), and that as t age increases, Cy (T) evolves towards 
the standard glass value (Fig. [I}. Inspection of the curves 
for large t age (Fig[2fc)) shows that a small pre-peak ap- 
pears at T « 113 K. If the standard protocol (Fig. 1) 
for the identification of the glass transition in the spe- 
cific heat is applied, a T g of » 75 K is derived (Fig|2fc)). 
The amplitude of the peak in Cy(T) is of the order of 
1 — 2.5 J/mol/K, and is reminiscent of the experimental 
value obtained in DSC measurements of hyperquenched 
water after annealing. In the present case, in which crys- 
tallization does not interfere with the heating scan, there 
is no ambiguity in associating this peak with a precursor 
of the true glass transition which takes place at a much 
higher T. To prove that the weak Cy pre-peak is outside 
the noise level, we show in Fig. [21 the T-dependence of 
the energy, where a clear maximum can be observed. 

Results presented in Figs. |21 and are consistent 
with recent anneal-and-scan experiments |l2j on hyper- 
quenched inorganic glass which does not crystallize on 



heating. Our simulations thus suggest that the mea- 
sured specific heat peak (Ref. |ll|) which has been used 
to identify T g is a pre-peak associated with the use of a 
hyperquenched sample combined with the annealing pro- 
cedure. 

Next, we discuss the possibility of mod eling th e simu- 
lation results using the TNM approach [lix l22t l25| . which 
is able to model the experimental heating scan of the 
specific heat for glasses generated with standard cooling 
rates (although it fails when applied to hyperquenched 
glasses 23] ). The TNM model assumes the response func- 
tion of the system can be represented by a stretched ex- 
ponential function with stretching parameter (3. It also 
assumes that the relaxation time r depends not only on 
the bath temperature T but also on a fictive temperature 
Tf which accounts for the out-of-equilibrium condition. 
Narayanaswamy proposed that r is related to Tf by 



r(T, Tf) = Aexp 



xAh* 
RT 



{l-x)Ah* 
RTf 



(1) 



where < x < 1, Ah* and A are constants and R is the 
ideal gas constant 26]. In fact, Eq. QJ corresponds to 
Narayanaswamy's original equation rewritten by Moyni- 
han [15| who introduced the parameter x. An alterna- 
tive relation^ is offered by the generalization of the 
Adam-Gibbs expression, which connects, in equilibrium, 
r to the configurational entropy S c . The resulting Adam- 
Gibbs-Scherer (AGS) expression is 



t(T, Tf) = A'exp 



E A 



Sc(T f )T 



where 



Sc(T) 



AC 
~T~ 



dT, 



(2) 



(3) 



Tr- is the Kauzmann temperature, Ea is a constant, and 
AC is the difference between the specific heats of the 
liquid and the glass. The TNM model requires (3 as a 
fitting parameter; additionally it requires the parameters 
(A, x, Ah*) for the Narayanaswamy-expression, or (A', 
Ea, Tk) for the AGS expression. 

Figure 21 compares our MD results during the heat- 
ing scan of the standard and hyperquenched glass with 
the predictions of the TNM model using both the 
Narayanaswamy and AGS expressions for t. A detailed 
analysis will be presented elsewhere 0; here we show 
that both the Narayanaswamy and AGS expressions give 
satisfactory results only for the standard glass (Fig.|3Ja)). 
The quality of the fit for the hyperquenched glass is un- 
satisfactory, as observed in the analysis of experimental 
data for hyperquenched samples 23], suggesting that in 
the hyperquenched experiments, the aging sample can- 
not be connected to a liquid at a fictive temperature 
Tf 29] . It also suggests that the application of the TNM 
approach for testing the shadow glass transition must be 
taken with caution|ll|. 



3 



In summary, we show by a proper numerical proto- 
col that the complex phenomenology of the glass transi- 
tion can be reproduced in simulations, notwithstanding 
the large differences in cooling rates. The TNM model 
which is able to describe the experimental specific heat 
for the standard cooling rate also describes the corre- 
sponding simulation results. One advantage of simula- 
tions is to shed light on phenomena that occur outside 
the experimentally-accessible region[30j. Our simulations 
show that the glass transition is characterized by a large 
specific heat peak, and when hyperquenched samples are 
annealed, the glass peak is anticipated at lower T by 
a pre-peak, with a much lower amplitude. This result 



supports the recent reinterpretation of the existing ex- 
perimental data , which identifies 136 K as the tem- 
perature of the pre-peak, and suggests « 165 K as the 
'true' glass transition temperature. 
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FIG. 1: Specific heat from MD simulations calculated by dif- 
ferentiating the total energy during heating of the standard 
glass (SG). Circles denote equilibrium values of CV(T) in the 
liquid state. 
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FIG. 2: (a) T-dependence of Cv for heating scans of the 
standard glass (SG) and the hyperquenched glass (HG). Also 
shown are the heating scans of the HG which has been an- 
nealed at T age = 100 K for four different aging times t age . 
(b) Magnification of (a). Curves are shifted for clarity by 
3 J/mol/K (for t age = 20 ps), 5 J/mol/K (for t age = 300 ps), 
7 J/mol/K (for t age — 1 ns), 9 J/mol/K (for t age — 20 ns), and 
11 J/mol/K (for SG). (c) Magnification of C V (T) for the hy- 
perquenched glass annealed at T age = 100 K for t age — 20 ns, 
to highlight the weak pre-peak at T f» 113 K. The straight 
lines show a possible construction that, in the absence of the 
peak at T ps 220 K, could be interpreted as the glass transi- 
tion temperature "T g " . 
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FIG. 3: T dependence of AE = E - 6RT, where E is the 
total energy per molecule, R is the gas constant, and 6RT is 
the total energy of a glass of rigid molecules in the harmonic 
approximation, contributing a constant 6R to CV ■ Using AE 
instead of E amplifies the very weak signal, whose derivative 
is responsible for the weak peak in CV shown in Figl^c). 
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FIG. 4: Heating scans from the (a) standard glass (SG) 
and (b) hyperquenched glass (HG) . Our simulations are com- 
pared with the predictions of the TNM model using both 
the Narayanaswamy model (fitting parameters for the SG are 
]n(A/ns) = -22.36, (3 = 0.525, x = 0.635, and Ah*/R = 
4632 K) and the Adam-Gibbs-Scherer expressions (fitting pa- 
rameters for the SG are \n(A/ns) = —9.86, (3 = 0.519, and 
Ea = 27626 kj/mol; we use S c from Ref. [32] and hence we 
do not require Tk)- We see that the TNM model describes 
the behavior of CV(T) for the SG but fails for the HG. 



